#include "uvp.h"



void calculate_dt(double Re,
                  double tau,
		  double *dt,
                  double dx,
		  double dy,
		  int imax,
		  int jmax,
		  double **U,
		  double **V)

{

	double a = (Re / 2.0) * (1.0 / ( 1.0 / pow(dx, 2.0) + 1.0 / pow(dy,2.0) ) ) ;
	double u_max = 0.0;
	double v_max = 0.0;
	int i , j;
	double b = 0.0 , c = 0.0;

	/* handle negative tau */
	if ( tau <= 0.0 )
		return;

	for ( i = 0; i < (imax + 1); i++) {
		for ( j = 0; j < (jmax + 1); j++) {
			u_max = fmax(u_max, fabs(U[i][j]));
			v_max = fmax(v_max, fabs(V[i][j]));
		}
	}

	/* the following checks are needed to avoid division-by-zero */
	if (u_max != 0.0)
		b = (dx / u_max);

	if (v_max != 0.0 )
		c = (dy / v_max) ;

	/* Now we need to check if the minimum value was found in U & V or b and c still have the initial values which are 0.0 */

	*dt = a;

	if ( b != 0.0 ) {
		*dt = fmin (*dt , b);
	}
	if ( c != 0.0 ) {
		*dt = fmin (*dt , c);
	}

	*dt = (tau) * (*dt);
}


void calculate_fg(double Re,
		  double GX,
		  double GY,
		  double alpha,
		  double dt,
		  double dx,
		  double dy,
		  int imax,
		  int jmax,
		  double **U,
		  double **V,
		  double **F,
		  double **G,
		  int **Flag
		  )
{
	int i,j;

	/* initialize the boundary values */
	for (j = 1; j <= jmax; j++) {
		F[0][j] = U[0][j];
		F[imax][j] = U[imax][j];
	}
	for ( i = 1; i <= imax; i++) {
		G[i][0] = V[i][0];
		G[i][jmax] = V[i][jmax];
	}


	/* calculate the interior values of F & G */
	for (i = 1; i < imax; i++) {
		for (j = 1; j <= jmax; j++) {
			if (is_fluid(Flag[i][j]) && is_fluid(Flag[i + 1][j])) {
				F[i][j] = U [i][j] + dt * (diffusion(Re,i,j,dx,dy,U) - \
					calculate_convection_x(i,j,dx,dy,alpha,U,V) + GX);
			}
		}
	}
	for (i = 1; i <= imax; i++) {
		for (j = 1; j < jmax; j++) {
			if (is_fluid(Flag[i][j]) && is_fluid(Flag[i][j + 1])) {
				G[i][j] = V [i][j] + dt * (diffusion(Re,i,j,dx,dy,V) - \
					calculate_convection_y(i,j,dx,dy,alpha,U,V) + GY);
			}
		}
	}

	/* Handle F & G on the obstacle boundaries */
    	for ( i = 1; i <= imax; i++ ) {
		for ( j = 1; j <= jmax; j++ ) {
			switch (Flag[i][j]) {
				case B_N:
					G[i][j] = V[i][j];
					break;
				case B_W:
					F[i - 1][j] = U[i - 1][j];
					break;
				case B_S:
					G[i][j - 1] = V[i][j - 1];
					break;
				case B_E:
					F[i][j] = U[i][j];
					break;
				case B_SW:
					F[i - 1][j] = U[i - 1][j];
					G[i][j - 1] = V[i][j - 1];
					break;
				case B_SE:
					F[i][j] = U[i][j];
					G[i][j - 1] = V[i][j - 1];
					break;
				case B_NE:
					F[i][j] = U[i][j];
					G[i][j] = V[i][j];
					break;
				case B_NW:
					F[i - 1][j] = U[i - 1][j];
					G[i][j] = V[i][j];
					break;
				default:
					break;
			}
		}
	}

}

void calculate_uv(double dt,
		  double dx,
		  double dy,
		  int imax,
		  int jmax,
		  double **U,
		  double **V,
		  double **F,
		  double **G,
		  double **P,
		  int **Flag
		  )
{
	int i,j;


	for (j = 1; j <= jmax; j++){
		for (i = 1; i < imax; i++){
			if (is_fluid(Flag[i][j]) && is_fluid(Flag[i + 1][j])) {
			/*if ( is_interior(Flag[i][j]) ) {*/
				U[i][j] = F[i][j] - (dt/dx)*(P[i + 1][j] - P[i][j]);
			}
		}
	}

	for (j = 1; j < jmax; j++){
		for (i = 1; i <= imax; i++){
			if (is_fluid(Flag[i][j]) && is_fluid(Flag[i][j + 1])) {
			/*if ( is_interior(Flag[i][j]) ) { */
				V[i][j] = G[i][j] - (dt/dy)*(P[i][j + 1] - P[i][j]);
			}
		}
	}

}


void calculate_rs(
		  double dt,
		  double dx,
		  double dy,
		  int imax,
		  int jmax,
		  double **F,
		  double **G,
		  double **RS,
		  int **Flag
		 )
{
	int i,j;
	for ( j = 1; j <= jmax; j++){
		for (i = 1; i <= imax; i++){
			if (is_fluid(Flag[i][j])) {
				RS[i][j] = (1.0/dt)*((F[i][j] - F[i - 1][j])/dx + (G[i][j] - G[i][j - 1])/dy);
			}
		}
	}
}


/*****************************************************************************************

The following equations are used in calculate_fg but were not specified in the worksheet.

This part is full of math. Let's hope that we get it right from the 1st shot!

*****************************************************************************************/

inline double calculate_convection_x(int i,
			     int j,
			     double dx,
			     double dy,
			     double alpha,
			     double **U,
			     double **V
			    )
{
  double c1,c2; /* c1 represents the convection of x in x-direction and c2 represents the convection of x in y-direction */
  c1 = (1.0/dx)*(pow((U[i][j] + U[i + 1][j])/2.0, 2.0) - pow((U[i - 1][j] + U[i][j])/2.0 , 2.0)) + \
       (alpha/dx)*(( \
       (fabs(U[i][j] + U[i + 1][j])/2.0) * ((U[i][j] - U[i + 1][j])/2.0)) - \
       (fabs(U[i - 1][j] + U[i][j])/2.0) * ((U[i - 1][j] - U[i][j])/2.0));

  c2 = (1.0/dy)*(((V[i][j] + V[i + 1][j])/2.0 )*((U[i][j] + U[i][j + 1])/2.0) - ((V[i][j - 1] + V[i + 1][j - 1])/2.0)*((U[i][j - 1] + U[i][j])/2.0)) + \
       (alpha/dy)*( \
		    (fabs(V[i][j] + V[i + 1][j])/2.0)*((U[i][j] - U[i][j + 1])/2.0) - \
		    (fabs(V[i][j - 1] + V[i + 1][j - 1])/2.0)*((U[i][j - 1] - U[i][j])/2.0)
		    );
  return (c1 + c2);
}

inline double calculate_convection_y(int i,
			      int j,
			      double dx,
			      double dy,
			      double alpha,
			      double **U,
			      double **V
			     )
{
	/* c1 represents the convection of y in x-direction and c2 represents the convection of y in y-direction */
	double c1,c2;
	c1 = (1.0/dy)*(pow((V[i][j] + V[i][j + 1])/2.0, 2.0) - pow((V[i][j - 1] + V[i][j])/2.0 , 2.0)) + \
       		(alpha/dy)*(( \
		(fabs(V[i][j] + U[i][j + 1])/2.0) * ((V[i][j] - V[i][j + 1])/2.0)) - \
		(fabs(V[i][j - 1] + V[i][j])/2.0) * ((V[i][j - 1] - V[i][j])/2.0));

	c2 = (1.0/dx)*(((U[i][j] + U[i][j + 1])/2.0 )*((V[i][j] + V[i + 1][j])/2.0) - ((U[i - 1][j] + U[i - 1][j + 1])/2.0)*((V[i - 1][j] + V[i][j])/2.0)) + \
	       (alpha/dy)*( \
		(fabs(U[i][j] + U[i][j + 1])/2.0)*((V[i][j] - V[i + 1][j])/2.0) - \
		(fabs(U[i - 1][j] + U[i - 1][j + 1])/2.0)*((V[i - 1][j] - V[i][j])/2.0)
		    );
	return (c1 + c2);
}


inline double diffusion(double Re,
		 int i,
		 int j,
		 double dx,
		 double dy,
		 double **M
		 )
{
	return (((M[i - 1][j] - 2*M[i][j] + M[i + 1][j])/(dx*dx)) + ((M[i][j - 1] - 2*M[i][j] + M[i][j + 1])/(dy*dy)))/Re;
}


